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ABSTRACT 

The structures produced during the epoch of reionization by the action of radiation 
on neutral hydrogen are in principle different to those that arise through gravitational 
growth of initially small perturbations. We explore the difference between the two 
mechanisms using high resolution cosmological radiative transfer. Our computations 
use a Monte Carlo code which raytraces directly through SPH kernels without a grid, 
preserving the high spatial resolution of the underlying hydrodynamic simulation. Be- 
cause the properties of the first sources of radiation are uncertain, we simulate a range 
of models with different source properties and recombination physics. We examine 
the morphology of the neutral hydrogren distribution and the reionization history in 
these models. We find that at fixed mean neutral fraction, structures are visually most 
affected by the existence of a lower limit in source luminosity, then by galaxy mass 
to light ratio, and are minimally affected by changes in the recombination rate and 
amplitude of mass fluctuations. We concentrate on the autocorrelation function of the 
neutral hydrogen, S^Hiij") as a basic quantitive measure of Radiation Induced Structure 
(RIS.) All the models we test exhibit a characteristic behaviour, with S^hi becoming 
initially linearly antibiased with respect to the matter correlation function, reaching a 
minimum bias factor 6 ^ 0.5 when the universe is ~ 10 — 20% ionized. After this S^ui 
increases rapidly in amplitude, overtaking the matter correlation function. It keeps a 
power law shape, but flattens considerably, reaching an asymptotic logarithmic slope 
of 7 ~ —0.5. The growth rate of HI fluctuations is exponentially more rapid than 
gravitational growth over a brief interval of redshift Az 2 — 3. 
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1 INTRODUCTION 

In the standard cosmological model, the large-scale structure 
in the density field grows from small initial perturbations 
through the mechanism of gravitational instability. Statisti- 
cal measures of this structure can be used to both verify the 
growth mechanism (see e.g., Bernardeau et al. 2002 and ref- 
erences therein) and quantify the initial pertubations. A dif- 
ferent kind of growth of structure is expected when we con- 
sider the neutral hydrogen density field during the epoch of 
reionization (see e.g., the review by Loeb & Barkana 2001). 
In this case, bubbles of ionized material form first around 
bright sources and grow as the ionization fronts overlap un- 
til the universe is fully ionized. Statistical measures applied 
to this "Radiation Induced Structure" (hereafter RIS) can 
be used in a similar way to the gravitational instability pic- 
ture above, but this time to verify the process of reionization 
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and the nature of the sources of radiation. The effect of RIS 
is likely to be qualitatively different from that of gravity, 
and as a result the statistical signatures will be different. 
Our aim in this paper is to explore the differences, using 
ray traced simulations of reionization. We aim to both find 
out how reionization is different from gravity in the way 
it forms structure and how to use these differences to cat- 
egorize reionization scenarios. The statistical properties of 
RIS have been explored in many other works, e.g., using 
the power spectrum of HI fluctuations by Furlanetto et al. 
(2004), Zaldarriaga et al. (2004), Morales & Hewitt (2004), 
and using Minkowski functionals by Gleser et al. (2006). In 
the present paper we will focus on the autocorrelation func- 
tion. 

Theoretical studies of reionization include both ana- 
lytic work, e.g., Miralda-Escude et al. (2000), Wyithe & 
Loeb (2003), Cen (2003), Liu et al. (2004), Furlanetto et al. 
(2004), and numerical simulations, e.g., Razoumov & Scott 
(1999), Abel et al. (1999), Gnedin (2000), Ciardi et al. 
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(2001), Sokasian ot al. (2001), Razoumov et al. (2002), 
Sokasian et al. (2004). 

Recently, N-body simulations with radiative transfer 
post processing have been performed by Iliev et al. (2006), 
and Zahn et al. (2007) with box sizes as large as 100 /t^^Mpc 
and able to resolve halos down to masses of ~ 2 X 10^ Mq. 
Kohler et al. (2007) have performed simulations with an ex- 
tremely large box size (up to 1280 /i~^Mpc) in which hydro- 
dynamics and radiative transfer are coupled self consistently. 
These simulation rely on higher resolution simulations to 
calibrate the sub grid (< 10 /i~^Mpc) physics. Other re- 
cent advances in simulating reionization include the work of 
Trac & Cen (2006), a hybrid N-body dark matter and RT 
approach, a code which was used to study the growth of 
bubbles during reionization by Shin et al. (2007). 

Analytic work has been carried out using perturbation 
theory to predict inhomogeneities in the density of neutral 
hydrogen and photons (Zhang et al. 2007), as well as Press- 
Sehechter based analyses (e.g., Furlanetto et al. 2004.) Most 
relevant to the work here, the autocorrelation function of 
21cm emission has been examined in analytic models by 
Wyithe & Morales (2007) and Barkana (2007). 

In the present paper we do not aim to simulate partic- 
ular observational probes of this epoch, such as the Doppler 
scattering of CMB photons on relativistic electrons, or the 
21cm emission from neutral hydrogen. Instead, we concen- 
trate on the differences between ^(r) for the density field 
and the density field modulated by RIS. In principle, this 
will be directly observable in the future, through the vari- 
ous observational probes (e.g., Caxilli et al. 2004, Peterson 
et al. 2005, Valdes et al. 2006) We leave exploration of this 
to future work. 

The format of this paper is as follows. In §2 we describe 
the simulations, including the N-body outputs and the ra- 
diative transfer code. We also describe our different models 
for the sources of ionizing radiation and different physical 
conditions we simulate. In §3 we show how global proper- 
ties such as the mean ionized fraction evolve as a function 
of redshift in the different runs, and then in §4 examine 
the morphology of the neutral and ionized hydrogen density 
field. In §5 we measure the autocorrelation function in our 
different models and explore how it can differentiate between 
them. We summarise our results and discuss them in §6 



2 SIMULATIONS 

We work in the context of the standard cosmological con- 
stant dominated universe, with parameters Qa = 0.7, 
flm = 0.3 fib = 0.04, and a Hubble constant Ho = 
70 km s~^ Mpc~^. The initial linear power spectrum is 
cluster-normalised with a linearly extrapolated amplitude 
of as = 0.9 at ^; = 0. The radiative transfer is carried out as 
post-processing on N-body simulation outputs. 

2.1 N-body outputs 

We run our N-body simulations with the cosmological hy- 
drodynamic code Gadget (Springel et al. 2001). Our fidu- 
cial simulation is run in a 40 ft~^Mpc cubical volume. Iliev 
et al. (2006) have shown in tests using subvolumes of a larger 
(100 h~^Mpc) simulation that 30 /i~^Mpc is the smallest 



box side length for which the scatter in reionization histo- 
ries in different volumes is reasonably small. We use 256"^ 
dark matter and 256'^ gas particles. The mass resolution is 
therefore 6.05 x lO^'M© for gas particles and 3.93 x 10® M© 
for dark matter. We also run higher resolution models in 
tests, as detailed below. We do not include radiative cooling 
or star formation when computing the gas dynamics so that 
our simulations are similar to those carried out by Sokasian 
et al. (2001). We run the RT as postprocessing, so that there 
is no coupling between the hydrodynamics and radiation. In 
this respect, the gas serves to trace the dark matter distribu- 
tion closely (there is little difference at these high redshifts.) 
We chose our sources of radiation to be associated with dark 
matter haloes (see below). 

In addition to the fiducial run, we run another model 
with identical box size and particle number but with differ- 
ent random initial phases in order to roughly indicate the ef- 
fect of simulation cosmic variance. We also run a model with 
a different amplitude of mass fluctuations (0-8= 0.7.) For 
resolution tests, we run a sinmlation with 128^ gas and dark 
matter particles in a box of size 20 h'^Mpc ( the same mass 
resolution as the fiducial run) and a simulation with 256^ 
gas and dark matter particles in a box of size 20 ft-^Mpc 
( eight times better mass resolution than the fiducial run). 
All models wore started at z = 50 and run until z = 5.2. We 
output snapshots of the density field every 25 Myr, so that 
there are approximately 40 snapshot files per run. 

2.2 Radiative transfer 

After choosing models for the sources of radiation (see be- 
low) we carry out raytracing simulations of radiative transfer 
(RT) to study the evolution of the neutral hydrogen density. 
The code we use to do this carries out Monte Carlo RT to 
follow photon packets through the distribution of matter. 
The code is based on that used by Croft (2004) to study the 
fiuctuating radiation background field at lower redshift but 
incorporates time dependent RT in a Monte Carlo manner 
similar to the CRASH code of Maselli et al. (2003) (we actu- 
ally only treat hydrogen here, so we are in fact closer to the 
earlier work of Ciardi et al. (2001). The code in the present 
paper traces directly through the SPH particle kernels (see 
also Kessel-Deynet & Burkert, 2000, Susa, 2006, Yoshida 
et al. 2007, and Daqle et al. 2007, for non grid based radia- 
tive transfer, and Semelin et al. 2007, and Oxley & Woolf- 
son 2003 for SPH codes that trace through a Barnes-Hut 
tree) and so requires no regridding of the density field be- 
tween outputs. The spatial resolution of the RT is therefore 
in principle higher in dense regions than would be possible 
with a uniform grid. This approach is described in detail in 
Altay et al. (2007), where test problems arc carried out. Wc 
also outline some of the relevant features of the code briefly 
below. 

In the present paper wo model only the hydrogen com- 
ponent of the Universe (assumed to comprise 0.76 of the 
baryonic mass). We also do not explicitly follow the temper- 
ature evolution of the gas, beyond taking temperatures of 
particles to be lO'' K when they are ionized and Tcmb when 
they are not. Wc follow coUisional ionization, photoioniza- 
tion and recombinations using the rates given in Cen (1992). 
We randomly sample source photons from a power law dis- 
tribution, = v", of photon energies (more details on the 
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sources are given in section 2.3, below.) In the present work, 
photon packets are monochromatic and are emitted isotrop- 
ically from sources, again using a random number generator 
to pick directions. 

Every time a ray is traced through a particle, the num- 
ber of recombination photons which have been produced in 
that particle since the last time it was visited are added to 
a stack. When the stack size reaches one packet, a recom- 
bination photon packet is emitted from the particle where 
this occurs. The frequency of the recombination radiation is 
given by the Milne relation (Osterbrock, 1989.) 

We have one numerical code parameter, c/ which sets 
the size of the photon packet, and therefore the time reso- 
lution of the code. This parameter c/ is the number of fully 
neutral simulation gas particles which could be ionized by 
one photon packet if the energy in that packet was split up 
into 13.6 eV photons. For example, if a packet consists of 

13.6 eV photons and a gas particle contains Nhi neu- 
tral hydrogen atoms then c/ = N^/Nhi- By trying runs 
with different values of c/ we change how well the code can 
resolve the recombination timescale by governing the av- 
erage interval between rays visiting particles and updating 
their ionization states. Packets of recombination radiation 
as well as source radiation are governed by c/ so that shot 
noise arising from the discreteness of recombination mod- 
elling can be controlled. We choose a suitable value of c/ by 
carrying out convergence tests (see e.g., §5.1). In practice we 
find that cj = 0.33 is adequate for modelling clustering of 
HI and is the value we use in our fiducial simulations. This 
results in from 0.5 — 1.0 x 10* photon packets being used in 
each of the runs. 

The paper Altay et al. (2007) is mainly concerned with 
presenting a closely related code, the publically available 
code SPHRAY, which has many additional features to the 
one used in this paper, including the ability to model tem- 
perature evolution and Helium species. 

2.3 Sources/runs 

Eventually, observations of the RIS at the epoch of reioniza- 
tion will be useful as probes of the sources of radiation as 
well as cosmology. The statistical measures of this structure 
are likely to be correlated with the properties of the sources, 
their luminosities, lifetimes, and clustering. One of the goals 
of this paper is to study the RIS produced by various ex- 
tremely simple models for the ionizing source population, 
in order to see how they can be differentiated (principally 
through the autocorrelation function, which we focus on) 
and which features appear to be generic to the models con- 
sidered. 

The range of possible sources for reionization is ex- 
tremely wide, including decaying dark matter (e.g., Mapelli 
et al. 2006), primordial black holes (Ping & Fang 2002, Ri- 
cotti et al. 2007), high redshift miniquasars (e.g., Madau 
et al. 2004), population III stars (e.g., Sokasian et al. 2004), 
population II stars (e.g., Sokasian et al. 2003), some more 
or less likely than others. Rather than attempting to sim- 
ulate particular models in detail, we restrict ourselves to 
simply parametrized models which relate the ionizing ra- 
diation intensity directly to the dark matter distribution. 
This is on the understanding that in most reasonable mod- 
els of reionization there would be some relationship between 



the two (either through galaxies and stars associated with 
dark matter overdensities, or directly through dark matter 
clumps decaying to ionizing photons). In particular, we as- 
sociate sources of radiation to dark matter halos. 

This approach has been used also by Mellema et al. 
(2006), who use a constant mass to light ratio to assign ion- 
izing radiation to each halo, as well by Zahn et al. (2007), 
who populate each halo with a single ionizing source whose 
luminosity is proportional to host halo mass. McQuinn et al. 
(2006) have also simulated 17 different variations of this type 
of model with various ionizing photon efficiencies, prescrip- 
tions for feedback and and minihalos. McQuinn et al. (2006) 
focus on the morphology of HII regions. 

All the runs we use for the main studies in this pa- 
per have a simulation box length of 40 /i~^Mpc and parti- 
cle number of 2 x 256^, although as explained in §2.1, for 
resolution studies we have some runs with different mass 
resolutions and box sizes. We find halos using a standard 
friends-of-friends routine, with a linking length of 0.2 times 
the mean interparticle separation. The minimum halo mass 
we use as a source in our fiducial run is 1.6 X lO'^h-iMo, 
containing only 8 (gas and dark matter) particles. This is ap- 
proaching the scale of mini halos expected to host numerous 
weak ionizing sources and provide small scale dumpiness to 
the IGM, but as with the calculations of Zahn et al. (2007) 
(who do have a smaller particle mass) it is still approxi- 
mately an order of magnitude too large. Although we only 
resolve such halos with a small number of particles, we have 
checked using simulations with 8 times better mass resolu- 
tion (see §5.1) that this does not affect our computation of 
the autocorrelation function of neutral hydrogen, the statis- 
tic we focus on. In general, the limited mass resolution of 
simulations will affect results both through the absence of 
sinks (minihalos) and sources hosted by small halos. Recent 
simulations (e.g., Santos et al. 2007) are beginning to ad- 
dress this directly through increases in simulation particle 
number. We return to this point in §5.1 and §6.2. 

We use the same density field for 10 of the 12 main runs, 
carrying out the RT as postprocessing using different source 
prescriptions. The other two are a low ffuctuation amplitude 
model (erg = 0.7) run with the same random phases and 
another model with the same amplitude as the fiducial case 
but with different phases. The 12 runs are differentiated by 
their different halo mass to light ratios, treatment of the 
relationship between halo mass and ionizing radiation, the 
recombination rate, and spectrum of radiation. An overview 
is given in Table [T] and they are descibed in detail below. 
We label the different simulations by short descriptive names 
rather than numbers or letters in order to avoid the necessity 
of the reader referring back to a table when examining the 
results. 

• The fiducial run. Here the instantaneous luminosity 
of sources is proportional to the dark matter halo mass, 
with L = Lo X Mhaio/h~^M0 In our fiducial run, we 
take I/O = 2.7 x lO"^^erg/s/h~^M0, and a source spec- 
trum with Fi, (X I'"*, appropriate for population II stellar 
sources (e.g., Sokasian et al. 2003). This simple model is also 
used by Zahn et al. (2007), who assume a number of pho- 
tons proportional to halo mass, (with a conversion factor 
3.1 X 10'*^ photons/sec/h~^MQ .) This corresponds to our 
fiducial model having a luminosity per halo mass 4 times 
larger than that of Zahn et al. . Using the computation of 
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Table 1. Radiative transfer simulations discussed in this paper. Further details are given in section 2.3 



Simulation (Luiniuosity Spectrum Recomb- Box length Particle Comments 

integrated to a -ination rate ( /i~^Mpc) number 

2 = 6) /fiducial (F„ oc i/") 



fiducial 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 




L/2 


0.5 


-4.0 


1.0 


40.0 


2 X 256^ 




L/4 


0.25 


-4.0 


1.0 


40.0 


2 X 256-^ 




L/8 


0.125 


-4.0 


1.0 


40.0 


2 X 256^ 




L indcp M 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 


luminosity independent of halo mass 


M > 10l°h-lMo 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 


lower limit on halo mass for sources 


M < 10l°h-lMo 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 


upper limit on halo mass for sources 


as = 0.7 


0.87 


-4.0 


1.0 


40.0 


2 X 2563 


low amplitude of mass fluctuations 


no recomb. 


1.0 


-4.0 


0.0 


40.0 


2 X 256^ 




2 X recomb 


1.0 


-4.0 


2.0 


40.0 


2 X 2563 




spectrum 


1.0 


-2.0 


1.0 


40.0 


2 X 256^ 




other ran. phases 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 


same params. as flducial, different ICs 


hires 


1.0 


-4.0 


1.0 


20.0 


2 X 256^ 




small 


1.0 


-4.0 


1.0 


20.0 


2 X 128^ 




different c/ 


1.0 


-4.0 


1.0 


40.0 


2 X 256^ 


3 runs: photon packet c/ = 3.0, 1.0,0.1 



Zahn et al. , the source output of our fiducial model can be 
considered to be very roughly equivalent to one with Popu- 
lation II stars forming with an efficiency of /, =0.1 from a 
Salpotcr IMF, with a stellar lifetime of At = 5 x IC yrs and 
an escape fraction of /esc = 0.04. The recombination rate in 
the fiducial simulation is taken to be that computed directly 
from the gas density, and diffuse recombination photons are 
treated. 

• Runs L/2, L/4 and L/8. These are the same as the 
fiducial run in every respect except that the ionizing lumi- 
nosity has been reduced by an overall factor of 2,4 or 8. 
These models can be considered to be equivalent to for ex- 
ample reducing the efficiency of star formation, and/or the 
ionizing escape fraction with respect to the fiducial run. The 
L/4 model therefore corresponds to the model simulated by 
Zahn et al. (2007). 

• L indep M Wc use the same halo list as the fiducial 
run, but instead of a luminosity proportional to halo mass, 
we assign the same luminosity to all halos. The total inte- 
grated luminosity to ^; = 6 is set to be equal to that in the 
fiducial run. 

• M > lO^*' h-^Mo, M < 10^" h^^Mo For these two 
runs, the same halo source list is used as in the fiducial run, 
but with either an upper or lower cutoff applied in the mass 
of a halo which can host a source. As in the previous model 
above, the total integrated luminosity to z = 6 is set to 
be equal to the fiducial run. These runs can be considered 
to roughly model the effects of feedback which might cause 
disruption of galaxy sources in small halos. 

• as — 0.7 This run is the same as the fiducial run, 
except using a simulation which has a significantly lower 
amplitude of mass fiuctuations, leading to later halo forma- 
tion times. The source luminosity was proportional to the 
halo mass in the same fashion as for the fiducial run, also 
with L() = 2.7 X lO^^erg/s/h"^M0. The total integrated lu- 
minosity to « = 6 is therefore less than in the fiducial run. 

• No recomb, 2x recomb In these runs, the recombi- 
nation rate of ionzed hydrogen was either set to zero or dou- 



bled, and the number of recombination photon packets ad- 
justed accordingly. Otherwise, the runs are the same as the 
fiducial run. These models can be thought of as parametriz- 
ing the effects of changing the clumping factor of unresolved 
gas. 

• i/"^ spectrum A harder spectrum than the 

spectrum used for the fiducial run was used here. The to- 
tal number of ionizing photons integrated to z = 6 was 
set to be the same as the fiducial run, so that Lq = 
1.0 X 10^2 erg/s/h-^Mo. This run gives a rough indication of 
the effects of a harder spectrum than that of Population II 
stars, although it was not chosen to reproduce the spectrum 
of any particular source. For example, at the ZAMS, Pop III 
stars are expected to have a spectrum with ~ , with 
a ~ —1.3 close to the hydrogen ionizing edge (Tumlinson 
et al. 2003). For representative composite spectra of AGN 
(e.g., Telfer et al. 2002) a value of a = —1.8 is possible. The 
effect of the harder spectrum will be to allow the photons 
to penetrate further. The "preheating" effect of hard pho- 
tons is not treated here as we do not follow the temperature 
evolution of the gas explicitly. 

• Other random phases As a rough indicator of the 
effects of cosmic variance, this run has the same parameters 
as the fiducial run, but uses a different underlying density 
field realization. 

We evolve all the simulations to redshift z = 5.5, irre- 
spective of whether they have achieved full reionization by 
then. 



3 GLOBAL EVOLUTION 

Because we have many different simulation runs with differ- 
ent treatments of recombination and models for the sources, 
wc need to decide at which redsliifts to compare our mea- 
sures of clustering. The different runs reach 50% mean ion- 
ized fraction by mass, Xm at redshifts ranging from z = 9—6, 
so that comparing them at the same redshift will mean com- 
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Figure 1. Evolution of the mean mass weighted hydrogen neu- 
tral fraction (top panel) and ionized fraction (bottom panel) with 
redshift. Note that the bottom panel is on a log scale. We show 
results for the 12 different models desribed in §2.3. 

paring different stages of reionization. We will therefore con- 
centrate in the paper on snapshots taken at the same mean 
ionized fractions in the different runs, rather than at the 
same redshifts. 

In Figure [1] we show how the mean mass weighted neu- 
tral fraction (linear scale) and mean mass weighted ionized 
fraction (log scale) vary with redshift for the different runs. 
As expected, the run without recombinations reionizes first, 
reaching a 1% ionized fraction by redshift z — 13 and be- 
ing 1% neutral sd, z — 7.8. As reionization proceeds, the 
relationship between mean ionized fraction x,n and redshift 
is approximately exponential, Xm = e"'^"^'-*, (here Zi is 
the redshift when the model is fully ionized) for this and 
the other models. This is roughly true for all models, ex- 
cept for the model where only large galaxies (halo masses 
M > lO'^°h~'^M0) are sources, and the model with erg = 0.7, 
both of which have a substantially more rapid change in 
ionization with redshift, Xm = e~^'®'^~^'' . This can be 
explained by the fact that galaxies massive enough to be 
sources only form relatively late in these two models. 

The models where the luminosity per dark matter halo 
atom was varied, fiducial, L/2, L/4, L/8 reach the Xm = 0.5 
point later by Az = 1 for each halving of the luminosity. The 
slope of logXm vs z is very similar for Xm < 0.1 for these 
models but then changes as reionization proceeds. We note 
that the L/A model reaches Xm = 0.5 at z ~ 7, similar to the 
model of Zahn et al. (2007), which has a similar spectrum 
and source luminosity (see §2.3 above.) 

More physical insight can be gained by looking at two 
other quantities as a function of redshift, the number of 




z 



Figure 2. Evolution of number of recombinations per atom (top 
panel) and the photon mean free path (bottom panel) with red- 
shift. In the top panel, we use the filled stars to mark the points 
when the mean mass weighted ionized fraction Xm = 0.5 and the 
open stars to mark Xm = 0.99 (not all runs have been carried 
out to this point.) We show results for the 12 simulation runs 
described in §2.3. 



recombinations per atom and the photon mean free path, 
which are plotted in Figure (2] To compute the former, we 
divide the cumulative number of recombination photons by 
the total initial number of hydrogen atoms in the simulation 
volume. We plot a symbol on the curves at the point when 
the ionized fraction by mass reaches Xm ~ 0.5 and another 
at Xm = 0.99 (not all models reach this stage). Most of the 
models have below 0.5 recombination photons per atom by 
the time they have fully reionized, indicating that recombi- 
nations do not play the major role in the process of reion- 
ization. The curves for the different models tend to flatten 
off slightly for the last half of the reionization process. One 
striking feature of these curves is the importance of the am- 
plitude of mass fluctuations, as- As we would expect, the 
low amplitude model (erg — 0.7), having less clumping has 
a lower recombination rate. When compared at the point 
when Xm = 0.5, this model has experienced 3.2 times fewer 
recombinations than the fiducial model. This is greater than 
the ratio of erg for the two models (1.65), due to nonlinearity 
of clustering and the greater cosmic time in the low ampli- 
tude model to get to this point. 

The photon mean free path as a function of redshift 
is also shown in Figure [2] The mean path length between 
absorptions has been calculated by averaging over the most 
recent 5 x 10^ photon packets, a relatively small number. As 
a result, sawtooth features in the curves can be seen, which 
correspond to the times at which the new density field snap- 
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shots were inputted (every 25 Myr). We have tested with 
more widely spaced outputs (50 Myr) that this does not 
afTect the underlying curves. Another artifact of the mea- 
surement is the behaviour at large mean free path. When 
averaging the path length, we only include photon packets 
that have not wrapped more than once round the box. As a 
result, when the volume is close to fully ionized, the average 
is taken over a biased sample, those photons which have run 
into an absorber in less than a box length. The mean free 
path we plot therefore falls. 

The overall behaviour of the mfp is as we would expect, 
with a rapid rise for most models at early times, with a 
mfp corresponding roughly to the size of the HIT regions 
around sources. For example at redshift z = 14, when the 
ionized fraction is ~ 10~^, the mfp in the fiducial model 
is ~ 10fe~^kpc. In the spectrum case, the low cross 
section of neutral hydrogen for the hard photons means that 
the mfp starts ofi' nearly an order of magnitude larger, ~ 
300^~^kpc. The hardest photons will travel much farther 
than this mean, and there is a low level of ionization (ionized 
fraction Xm ~ 10"*) present throughout the whole volume, 
even at these very early redshifts. 

In order to understand how reionization proceeds, it is 
instructive to look at the relation between neutral fraction 
and density. This is plotted for the fiducial run in Figure [3] 
at 4 different stages in the evolution of the model. In each 
panel, we also show a histogram of particle density values 
and neutral fractions as well as a scatter plot of one versus 
the other. If we go from panel to panel, we can see that 
the high density regions are reionized first (the 'inside-out' 
scenario also found in Iliev et al. 2006). For example, when 
Xm = 0.1, the cloud of partially ionized particles to the 
right of the panel is centered around pj (p) ~ 10, whereas at 
Xm = 0.7 it is centered around pj (p) ~ 3. Within this cloud 
of points, the higher rate of recombinations in high density 
regions means that within the partially ionized volume there 
is a trend for higher neutral fractions at high densities. This 
trend is additional to the opposite trend for high density 
regions to be reach this partially ionized state first. As the 
universe comes close to being fully ionized, Xm = 0.97, par- 
ticles close to the mean density are the only ones still fully 
neutral and the particles with the highest densities have 
neutral fractions < 10~^. 

In the next section, we examine the morphology of the 
neutral and ionized regions in order to investigate this in 
more detail. For now, we can compare the scatter plots of 
neutral fraction and density in different models. We do this 
at Xm ~ 0.5 in Figure U The absence of a cloud of partially 
ionized points in the "no recomb" simulation is expected, as 
once particles are ionized, they drop off the bottom of the 
plot. Also as expected, the spectrum simulation and the 
"2x recomb" run have a greater density of particles in this 
region than the fiducial model. 



4 MORPHOLOGY 

Just as looking at structure in galaxy redshift surveys (e.g., 
Schectman et al. 1996) revealed filaments, voids and clusters, 
the RIS is expected to give rise to a complex morphology. 
In the case of structure in the density field, reproducing the 
visual characteristics of observational data was one of the 




-0.1 1 10 1000.1 1 10 100 



/>/</•> pJ </» 

Figure 3. Scatter plot of gas density (in units of the mean) 
against hydrogen neutral fraction for particles in the fiducial sim- 
ulation run (see §2.3). We show results for 4 different output 
times, characterized by the mean mass weighted ionized fraction 
Xm which appears in the panel labels. In each panel, we also show 
as shaded areas histograms of the number of particles in bins of 
hydrogen neutral fraction and also histograms binned by density, 
p/ < p >. The height of each histogram bin is on a log scale. 

drivers in searching for the correct theory of structure for- 
mation, and cosmological N-body simulations are expected 
to give rise to a "cosmic web" with the same appearance 
(e.g.. Bond et al. 1996). Deciding how to view the morphol- 
ogy of RIS is complicated by the question of whether to plot 
the neutral fraction, neutral density, ionized density or ion- 
ized fraction, each of which can lead to a potentially differ- 
ent impression. Also, unlike pure density fluctuations, which 
evolve relatively slowly, the exact time the RIS is plotted as 
reionization proceeds can yield very different results. Again, 
we will choose to compare different models at the same value 

of Xm- 

In Figure [5] we show a thin (1 /i~^Mpc) slice through 
the fiducial simulation volume, at the time when Xm = 0.5. 
We use a two dimensional color scale to show both the den- 
sity and the neutral fraction of hydrogen, as well as overplot- 
ting the positions of the individual sources of radiation, as- 
sociated with dark matter halos. It is apparent from the plot 
that low density regions where there are only a few isolated 
sources have not yet created noticeable bubbles, but that 
the highly clustered regions, associated with filaments and 
proto clusters have appreciable Stromgren spheres around 
them. Iliev et al. (2006) have measured the size of the bub- 
bles in their simulations by fitting spheres into the ionized 
regions, finding a median bubble radius of ~ 5 /i~^Mpc (see 
their Figure 13) at this late stage of reionization. Many other 
studies, e.g.. Shin et al. (2007), Zahn et al. (2007) have been 
carried out on the size of bubbles in simulations. Visually, 
our simulation appears to be broadly consistent with the 
sizes found by Iliev et al. , and we leave to future work de- 
tailed statistical characterization of the size and shape of 
voids in the neutral hydrogen. 

For now, we will comment on the obvious differences 
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Figure 5. A thin (1 /i~^Mpc) slice through the fiducial simulation volume (40 /i~^Mpc wide), at redshift z = 8.5, when the mean mass 
weighted H neutral fraction Xm = 0.5. We show both the density in units of the mean and the H neutral fraction using a two dimensional 
color scale. The positions of individual sources of ionizing radiation are also shown as (red) points. 



apparent between the morphology of RIS in the simulation 
and that which can be seen in the underlying density field. 
The RIS has a much higher contrast level, with the neutral 
fraction in HII regions being ~ 10^ times less than in the 
rest of the volume. The edges of HII regions are consequently 
much sharper than those of voids in the matter distribution, 
lending themselves to easier detection by void-finding tech- 
niques (e.g., Colberg et al. 2005). 

The matter density field over the redshift range relevant 
to reioniziation evolves little compared to the neutral density 
field. Reionization (say the change from Xm = 0.001 to Xm = 
0.999) in these models takes place over an change in scale 
factor a of ~ 2, and as £7™ ~ 1 to good approximation at 
these redshifts, linear growth of matter fluctuations occurs 
by the same factor. Because of this, most of the change in 
the morphology and structure of neutral density field occurs 
in the RIS. Plotting a slice through the neutral density field 
at different epochs allows us to see this well. In figure |6l 
we show this for 3 different models (Fiducial, M > 10^^, 
M < 10^"), at times when Xm varied from 0.1 to 0.9 in steps 
of 0.2. 

McQuinn et al. (2007) have shown that the large 
scale morphology of ionized hydrogen bubbles depends most 
strongly on Xm and the properties of the ionizing sources, 
and is relatively less affected by the specific subgrid model 
used to determine small scale source suppression and clump- 
ing factors. In Figure [6l we can see that models in the same 
column (same value of Xm) are indeed fairly similar. How- 



ever, because we use the mass-weighted ionization fraction 
Xm rather than the volume weighted fraction x^, our conclu- 
sions about the similarity of the morphologies is somewhat 
different than those of McQuinn et al. (2007). For exam- 
ple, the panel with Xm = 0.5 (middle panel) for the fiducial 
model (top row) seems to be most similar to the Xm = 0.3 
(2nd panel) for the M > 10^° model (middle row). This 
is because in the M > 10^° model, the ionized material 
is all concentrated in large bubbles, whereas in the fiducial 
model there are many smaller HII regions around the fainter 
sources which are hard to see but which nevertheless account 
for much of the mass in ionized material. 

A good way to see that this is the case is to refer to 
Figure [71 which shows the same slices through the same 
models, but this time plotting the ionized density rather 
than neutral density. The HII regions around the fainter 
sources are clearly seen in the top and bottom rows. The 
early HII regions appear sharper and perhaps more spherical 
when seen directly in terms of their ionized density, rather 
than as shadows in the neutral hydrogen plot (Figure [6]). 
This is understandable because of the fact that the edges 
of the bubbles and the totality of bubbles that are smaller 
than the slice thickness (1 h~^Mpc) in size will be somewhat 
obscured in Figure [6] by neutral hydrogen that lies in front 
or behind the bubble but is still in the slice. 

In Figures [7] and |6] we see little difference between the 
M < lO^^h'^MQ model and the fiducial model, indicating 
that the absence of the most massive halos does not greatly 
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Figure 6. Slices (1 h^^Mpc thick) through the neutral H density field in 3 simulation runs as a function of mean mass weighted ionized 
fraction, Xm- The top row shows the fiducial run, the middle row the M > lO^''h~^M0 run and the bottom row the M < lO^^h^^M© 
run (see §2.3 for full descriptions.) In each panel we use a log scale, with light shades representing high neutral H density. From left to 
right, we show results for Xm = 0.1,0.3,0.5,0.7 and 0.9. 



afTect the morphology, as long as we compare at the same 
value of Xm- If we look at the panels at the farthest right of 
these plots, the end stages of reionization [xm = 0.9), we can 
see that the models with small halos allowed do have HI rem- 
nants with more ragged edges and more small scale structure 
apparent in them than the M > lO^°h~^M0 model. In fu- 
ture work, it would be interesting to investigate the mass 
function of the disconnected HI remnants present at these 
times as they may have constraining power (as well as being 
likely detectable in 21cm emission. The Xm value at which 
the ionized regions percolate seems likely to depend on the 
source model also, as for example the Xm = 0.5 panel for the 
M > lO^^h'^Mo model consists less of disconnected HII re- 
gions than the other two panels. This makes sense, as the 
ionized material is more clustered, being closer to the more 
massive haloes. 

This Xm = 0.5 epoch is therefore a good one at which to 
compare the morphology of the other models also. In Figure 
[S] we show the same thin slice through the neutral density 
field for the 12 different models. The panel most different 
from the others is not suprisingly that for the "other ran- 
dom phases" model, for which the slice intersects a fairly 
spherical bubble (at the bottom of the plot) and a large 
region which has little sign of reionization (in the middle). 
Comparing this panel to the fiducial model, it might seem as 
though the latter is closer to the percolation stage, although 
they both have the same Xm = 0.5. In the later stages of 
reionization, when the ionized structures are a substantial 



fraction of the box size, we must naturally be careful in our 
interpretation of the morphology, due to the effects of cosmic 
variance and finite volume. As stressed by Iliev et al. (2006) , 
and others, large boxes are necessary to capture these pro- 
cesses, particularly at the later epochs. 

Looking systematically at the different panels of Fig- 
ure [8l we can see that in the fiducial, L/2, L/4 and L/8 
models as the luminosity per halo atom is reduced, the vi- 
sual extent of the structures becomes smaller and smaller. 
From Figure [2] we can see that in these models, the num- 
ber of recombinations becomes progressively higher as we 
reduce the luminosity, and although one would expect that 
this would be compensated for in the morphology by the 
fact that we compare all models at the same Xm, there is 
in fact still a large difference. We can therefore state that 
changing the luminosity of sources at fixed Xm does affect 
the visual morphology quite strongly. For the L/8 model, 
much of the ionized hydrogen is again hidden from view in 
the strongly neutral regions, as happened in the comparison 
of the fiducial and M > lO^°h~^M0 models. From a rough 
visual comparison of the different columns of Figure [6] and 
models L/2-L/8 in Figure [8] it seems as though scaling the 
luminosity of sources by a factor of ~ 4 changes the visual 
morphology of the HI density field in a similar fashion to 
changing Xm by ~ 0.2. Of course this will not hold for the 
morphology of the ionized regions, and we will investigate 
this statistically when we look at the correlation function in 
the next section. 
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Figure 7. This figure is the same as Fig. [6] except that we show the ionized density: plotted are shces (1 /i~^Mpc thick) through the 
ionized H density field in 3 simulation runs as a function of mean mass weighted ionized fraction, Xm- The top row shows the fiducial run, 
the middle row the M > 10^'^h~^MQ run and the bottom row the M < lO^^h^^MQ run (see §2.3 for full descriptions.) In each panel 
we use a log scale, with light shades representing high ionized H density. From left to right, we show results for Xm = 0.1,0.3,0.5,0.7 
and 0.9. 



5 AUTOCORRELATION FUNCTION 

Quantitative comparisons of the structure in the various 
models can be carried out by looking at the autocorrelation 
function, £^r- We will focus on ^(r) measured for the neutral 
density distribution, although we will briefly examine ^(r) 
for the ionized density field. We note that the Fourier trans- 
form of ^{r), the power spectrum P{k) of fluctuations has 
been studied in reionization models by many authors. We 
compute ^(r) for the gas density directly from the particle 
positions in the simulation, and for ^Hi{r) we weight each 
particle by its neutral fraction: 

^Hiir) = — -2 1, (1) 

iVp,e (XhI) 

where Np is the number of pairs of particles in a bin centered 
on r, Np^e is the expected number for a random distribution, 
{xHi)j and {xHi)k are the neutral fractions of particles in 
a pair, and {xhi) is the mean neutral fraction by mass. As 
with our examination of morphology of the neutral density, 
we will compute ^(r) for different simulation outputs chosen 
by their mass weighted ionized fraction, Xm- 



5.1 Resolution and boxsize tests 

In order to test the range of validity of our ^(r) results, we 
carry out several resolution and boxsize tests. Because the 



bubble like structures which overlap during the end stages of 
reionization occupy a large volume, one would expect that 
the simulation box size may have a strong effect on our re- 
sults. We therefore compute ^(r) for the neutral gas and the 
total gas density for two simulations with different box sizes 
(we vary the box side length by a factor of 2) but the same 
mass and spatial resolution (this is kept the same as our 
fidcucial model) . The results are show in Figure (5] where 
the 3 panels show €,{r)p and €,{r)Hi for at difi^erent stages of 
reionization (parametrized by Xm) in the fiducial model. 

We can see that the £,{r)p curves are extremely simi- 
lar over the range 0.05 /i~^Mpc < r < 4 /i~^Mpc for the 
3 values of Xm- In particular, a power law fit to the curves 
over this range gives virtually identical parameters. This is a 
very good thing, because it shows that no non-linear gravi- 
tational mode coupling has taken place with large scale den- 
sity modes of order the box size. This is one advantage of 
working at these high redshifts (z ~ 8 for Xm ~ 0.8 in this 
case), where we can see that even a 20 h~^Mpc box is large 
enough to study gravitational clustering over this range of 
length scales. The correlation function of the neutral gas 
density has a broadly similar behaviour. The position of the 
break in the powerlaw form of ^(r) is due to the finite size 
of the simulation volume and so we will restrict our analysis 
of ^(r) to smaller scales. 

Interestingly, the panel with results closest to the end of 
reionization {x,n) does not show any greater disagreement 
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Figure 8. Slices (1 h^'^Mpc thick) through the neutral H density field in the 12 simulation runs described in §2.3. We show all results 
at output times that correspond to a mean mass weighted ionized fraction, Xm = 0.5. The density is shown on a log scale (light colours 
for higher neutral density). Regions which are completely black generally have a neutral fraction < 10^^. 



for large r for (^{r)Hi than for £,{r)p. This means that at 
least over this limited range of scales, our measurements 
will also be reliable for the neutral density. On the smaller 
scales (r < 0.1 /i~^Mpc, there appears to be a rapid drop 
off in for the neutral density for both boxsizes. We shall 
see below when we consider mass and spatial resolution that 
our results will not be useful below these scales in any case. 

For our next test, we keep the box size fixed at 
20 /i~^Mpc but vary the particle number (and hence mass 



resolution) by a factor of 8. We also vary the spatial (force) 
resolution by a factor of 2. The coarser mass/spatial res- 
olution is the one we use in our fiducial case (only with 
a 40 /i~^Mpc box). The dark matter halos which we use to 
place our sources of radiation have the same lower mass cut- 
off 1.6 X 10^ h-^Mg in both runs, which corresponds to 8 
times fewer particles in the low resolution case. The results 
for ^(r)p and ^{r)Hi are shown in Figure [Till again for differ- 
ent values of Xm- The ^(r)p curves show good agreement on 
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Figure 4. Scatter plot of gas density (in units of the moan) 
against liydrogen neutral fraetion for particles in six different sim- 
ulation runs (taken from the 12 in §2.3). We show results for all 
simulations at the time when the mean mass weighted ionized 
fraction a;^ = 0.5 In each panel, we also show as shaded areas 
histograms of the number of particles in bins of hydrogen neutral 
fraction and also histograms binned by density, p/ < p >. The 
height of each histogram bin is on a log scale. 



large scales (the two simulations were run with the same ini- 
tial phases for the density field). On scales r < 0.15 /i~^Mpc, 
the two curves diverge, indicating the effects of mass and 
spatial resolution on the gravitational evolution of the gas 
density field. We note that on scales comparable to this the 
clustering in the gas density will be influenced by cooling 
and star formation, which we do not include here in any 
case. The £,{r)Hi correlation functions also agree well down 
to r < 0.15 /i~^Mpc, when they diverge even more sharply. 
We notice that as will the boxsize test, the worst disagree- 
ment on large scales occurs rather unexpectedly for the early 
stages of reionization [xm ~ 0.2). 

The flnal simulation parameter which we vary is c/, the 
maximum number of simulation particles which can be ion- 
ized by a single photon packet. This parameter (see Section 
2.2 for more details) is inversely proportional to the total 
number of photon packets used to carry out the Monte Carlo 
RT. For larger values of Cf, the radiation field will not be 
as smooth, and there will be more shot noise in the neutral 
density field. Our fiducial value of c/ = 0.3, and in Figure 
[11] we show what happens when this is varied from Cf — 3 
to Cf — 0.1, with all other simulation parameters the same 
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Figure 9. The effect of boxsize on ^: We show results for two 
simulations, the fiducial run (which has a 40 h~^Mpc cubical 
simulation volume, bold lines) and another with identical mass 
resolution and other parameters but in a 20 h~^Mpc box (thin 
lines.) Panels (a)-(c) show results for different outputs, labelled 
by the value of the mean mass weighted ionized fraction, Xm- As 
dashed and solid lines we show the autocorrelation function of 
the HI density and the gas density respectively. 
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Figure 10. The effect of mass and spatial resolution on We 
show results for two simulations, both in a 20 /i~^Mpc cubical 
volume. One has the same mass and spatial resolution as the fidu- 
cial run (meaning 128^ gas particles in this volume, thin lines) 
and the other has 8 times better mass resolution and two times 
better spatial resolution (meaning 256^ gas particles in this vol- 
ume, thick lines) Panels (a)-(c) show results for different outputs, 
labelled by the value of the mean mass weighted ionized frac- 
tion, Xm. As dashed and solid lines we show the autocorrelation 
function of the HI density and the gas density respectively. 
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Figure 11. The effect of the photon packet size parameter, Cf 
(sec §2.2) on ^. We show results for 4 simulations, the fiducial run 
(which has = 0.3) and 4 other runs with Cj^ = 3, 1 and 0.1. The 
fiducial run is shown as the thickest line and the others as thin 
lines. Smaller values of Cy imply a larger total number of photon 
packets. Panels (a)-(c) show results for different outputs, labelled 
by the value of the mean mass weighted ionized fraction, Xm- As 
dashed and solid lines we show the autocorrelation function of 
the HI density and the gas density respectively. 



as our fiducial model. The £,{r)p curves are almost identical, 
as might be expected, with the small differences due to the 
fact that reionization occurs at slightly different times in the 
different c/ runs. The curves for ^(r)^/ are also very simi- 
lar, and there is no apparent systematic effect even when c/ 
is made 10 times larger than our fiducial value. At least as 
far as the correlation function is concerned, we have there- 
fore converged to stable results with our fiducial number of 
photon packets. 

The tests in this section have therefore revealed that our 
results for S,p{r) should be reliable over at least the scales 
0.15 h~^yipc < r < 4 /i~^Mpc. We will concentrate on this 
range in our analysis, for example looking at the power law 
nature of £,p{r). For looking at larger scales, approaching 
the scale of bubbles at the time of percolation, larger sim- 
ulation volumes should be run in the future. We note that 
because £,p{r) for a 20 ft~^Mpc volume converged with a 
larger box on scales below r < 4 /i~^Mpc it is probably safe 
to assume that we can draw information from our fiducial 
volume (40 /i~^Mpc box) on scales up to r ~ 8 /i~^Mpc. 



5.2 The evolution of ^(r) 

We show ^(r) for the fiducial model in figure [12] for values 
of the ionized fractions Xm ranging from 0.1 to 0.99, which 
corresponds to a redshift range of 2 = 10.2 to z = 7.6. 
Before reionization starts, when Xm = 0, £,p{r) and £,Hi{r) 
are identical, by definition. However, once Xm has reached 
0.1, one can see that £,Hi{r) is somewhat lower than ip{r), 
by a factor of ~ 0.7, but has the same shape, a power law: 
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Figure 12. The autocorrelation function, ^(r) for the gas density 
field (solid lines) and the HI density field (dashed lines) in the 
fiducial simulation (see §2.3.) We show results for 11 different 
output times, for which we have labelled the HI curves with the 
mean mass weighted ionized fraction at that time, Xm- 



with slope 7 ~ 1.5 on scales r ^ 7 /i~^Mpc and a gentle 
break above it (at a scale dictated by the finite box vol- 
ume). We will explore power law fits to ^(r) in section 5.4 
below. For now, we explore qualitatively the behaviour of 
^p(r) as reionization proceeds. We note that this behaviour 
(a decrease and then an increase in the amplitude of clus- 
tering) has been seen in analytic calculations of the 21 cm 
brightness correlations (e.g., Wyithe & Morales, 2007) 

As Xm increases, ^p(r) exhibits the usual linear growth, 
with the amplitude increasing as expected under gravita- 
tional instability, and the shape remaning constant. As can 
be seen from the solid lines in Figure [121 this gravitational 
amplification of structure is very small (a factor ~ 1.5^) over 
the time of reionization. £,Hi{r) on the other hand shows 
dramatic growth, by a factor 100 or more over the same in- 
terval, as we would predict for example by examining the 
morphology of structure in the HI density field in Figure [6] 
The initial stage of reionization affects primarily the high 
density regions around sources of radiation. Removing their 
contribution from the clustering of HI leads to an antibias 
on all scales in (,Hi{r) with respect to S,p{r). After this, the 
effect of removing highly clustered neutral regions competes 
with the amplifying effect of RIS, with the latter winning 
after Xm ^ 0.35, raising ^Hi{r) above ^p{r) at this point. 
As the RIS grows in scale, the shape of ^Hi{r) begins to 
change, with the slope of the powerlaw region becoming flat- 
ter, reaching 7 ~ 0.5 when Xm = 0.99. 

The relationship between ^Hi{r) and ^p(r) can be ex- 
amined by plotting the bias as a function of scale, deflned 
by 



bHi{r) = 



^Hi{r) 



1/2 



(3) 



^(r) = (r/ro)" 



(2) 



. This quantity is plotted in Figure[T3]for the fiducial model, 
where it can be seen that bniij') is approximately constant 
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Figure 13. The HI bias factor, hni as a function of scale, r for 
the fiducial simulation. The factor bj://(r) is given by Equation |3] 



with scale for r ^ 5 h^^Mpc for Xm < 0.35, with a value 
less than 1 for low values of Xm- After this, bHi{r) takes 
on a positive power law slope. If ^///(r) and ^p(r) are fit by 
powerlaws, with slopes 7h/ and 7p, then we expect the slope 
of a power law fit to bni, 

bHiir) = {r/ro)-^'\ (4) 

to give 7f, = ^{■jHi — 7p). If we examine bni for Xm — 0.99, 
we find a slope 75 = i(0.5 — 1.5) — —0.5, as expected. 
We can see that this slope is reached quite rapidly after Xm 
reaches ~ 0.5 and does not appear to evolve much after this. 
We shall explore this further below. The behaviour of bni 
in this strongly ionized regime, is not constant with scale, 
indicating that during the late stages of reionization, 21cm 
tomography measurements (if they were possible on these 
small scales, which is observationally difhcult) would not re- 
sult in a straightforward measurement of the underlying gas 
or dark matter clustering. During the earlier stages, how- 
ever, the (anti)bias is linear, although its amplitude may be 
related to reionization process in a complex way. 

If instead we compute ^h+{i~)j the correlation function 
of the ionized hydrogren (e.g., as plotted in Figure . we 
find that the situation is simpler. This is shown for the fidu- 
cial model for different Xm values in Figure [T4l At the highest 
values of Xm, we obviously expect (r) ~ £,p(r), which has 
the characteristic power law behaviour for r ^ 8 h~^Mpc. At 
earlier times, ^£f+ (r) has the same slope, but a larger ampli- 
tude, indictating that the ionized regions trace a biased sub- 
set of the density distribution, around peaks in the density 
field (see e.g.. Kaiser, 1984.) When Xm = 0.1 (at z = 10.2), 
the value of ro is ~ 3 /i~^Mpc, (compare to 5 /i~^Mpc for 
L* galaxies at the present day, Zehavi et al. 2005). We can 
see from Figure [141 that Ch+('') breaks from a powerlaw at 
around r ~ 1 /i~^Mpc, which, from looking at Figure ?? 
corresponds approximately to the size of ionized bubbles at 
this early epoch. 
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Figure 14. Clustering of ionized hydrogen: we show the autocor- 
relation function, ^(r) for ionized H density field (dashed lines) 
and the gas density field (solid lines) in the fiducial simulation (see 
§2.3.) We show results for 11 different output times, for which we 
have labelled the ionized H curves with the mean mass weighted 
ionized fraction at that time, Xm. 

5.3 ^(r) for the different models of reionization 

We have seen that ^Hi{r) in the fiducial model has a similar 
shape to that which arises under gravitational instability. 
We now investigate (,Hi{r) in the other reionization models 
in order to see if there is a dependence on the properties 
of the sources of radiation or the physics of reionization. 
We again compute ^Hi{r) at output times corresponding to 
different specific values of Xm- 

The results for our 12 models are shown in Figure 1151 
along with power-law fits (described more fully in Section 5.4 
below) to the region 4.0 /i"^Mpc > r > 0.2 h'^Mpc. It is 
apparent that all models display behaviour broadly similar 
to the fiducial case, with a form reasonably approximate 
to a power law over ~ 1.5 decades in scale. As indicated 
by our tests above, the break on scales r ~ 7 /i~^Mpc is 
largely caused by the finite box volume and is expected to 
be similar in all cases. All models have a certain level of 
antibias between (,Hi{r) and (,p{r) at first, with the most 
extreme antibias being reached in the "L/8" model. The 
"M > lO'^°h~^M0" model has only a small level of antibias, 
indicating that the large bubbles formed rapidly by the very 
luminous sources quickly modulate the HI field. 

The correlation function becomes shallower and its am- 
plitude increases dramatically for all models towards the end 
of reioniziation. When Xm ~ 0.99, the models all have a very 
similar amplitude and a shallow power law slope 7 ~ 0.5. It 
should be noted that model "L/8" and "L indep M" do not 
have outputs for the very end stages of reionization. The lat- 
ter has very similar behaviour to the fiducial case, and the 
former is rather different, as we shall see when we consider 
the evolution of the power fit parameters. 

Overall one can see from the panels in Figure [TJlthat the 
RIS leads in all models to a strong increase in ^Hi{r) over 
the course of reionization. Unlike linear growth of fluctua- 



14 R.A.C. Croft and G. Altay 




1 10 1 10 1 10 

r (h"' Mpc^ r fh"' Mpc) r (h"' Mpc) 



Figure 15. The autocorrelation function ^ for all 12 models described in §2.3. The dashed lines are ^(r) for the neutral HI and solid 
lines for the gas density. The dotted lines show power law fits to the HI ^(r) for the region 4.0 /i~^Mpc > r > 0.2 h~^Mpc. For all panels 
except those for the "L/8 "and "L indep M" models, we show results for the same values of the mass weighted mean ionized fraction as 
in Figure ri2l i.e. Xm = 0.1 — 0.9 in steps of 0.1 and then Xm = 0.97 and Xm = 0.99. For the "L indep M" models, the Xm = 0.99 lines 
are not shown. 



tions under GI, the slope the correlations changes, reaching 
a similar shallow value in all models. No obvious features 
with any particular physical are created in £,Hi{r) in any of 
our variations of the reionization scenario. Despite the rel- 
atively wide variety of models employed, the differences in 
^Hi{r) are subtle, and will be explored in ro — 7 space below. 

5.4 Povirer law fits 

The galaxy-galaxy correlation function has only small devi- 
ations from power law behaviour over ~ 3 decades in length 
scale (e.g., Peebles, 1980, and Zehavi et al. 2002), despite 
the undoubted complexities of galaxy formation. The dark 



matter correlation function in simulations also exhibits this 
behaviour (e.g., Jenkins et al. 1998, Kravtsov et al. 2004.) 
It is therefore not unreasonable to expect that this type of 
scale invariance might also be created by RIS. We have seen 
from Figure [15] that this is indeed the case. In this subsec- 
tion we examine the power law behaviour of £,Hi{r) more 
quantitatively. 

We fit power laws to the region 4.0 /i~^Mpc > r > 
0.2 /i~^Mpc, based on the resolution and boxsize tests of 
Section 5.1. The power-law fits were carried out assum- 
ing Poisson errors on the (,Hi{r) points, so that the er- 
ror is (T = (1 -I- £,Hi{r)) / yJ\Np), where Np is the number 
of pairs of particles in a bin. We have tried changing A'p 
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Figure 16. Parameters for power law fits to the autocorrelation 
function of the HI density in all 12 models described in §2.3. Top 
panel: correlation length rg as a function of mass weighted ionized 
H density Xm- Middle panel: slope 7 as a function of Xm- Bottom 
panel: ro vs. 7. 



to include only pairs of particles above a threshold neu- 
tral fraction (e.g., xhi = 0.5), but this has a negible effect 
on the fit parameters. Also, changing the fitting region to 
1.0 h~^Mpc > r > 0.2 /i~^Mpc does not change any of our 
conclusions below. 

In the top 2 panels of Figure 1161 we show the depen- 
dence of the power law fit parameters ro and 7 on Xm for 
our different models. Looking at ro first, we can see that 
there is a strong dependence on Xm- The models all follow 
a trend roughly equivalent (within a factor of 2 in ro for 
all but one model) to logjQ ro ~ 0.06/(1. — Xm)^, with the 
M > 10i°h-^MQ model being the most extreme outlier. 
The fact that the lines for nearly all the different models 
are clustered together so tightly is rather surprising, and 
appears to be an indication that Xm plays a dominant role 
setting ro. At least in this fashion, the neutral fraction gov- 
erns the structure in the neutral hydrgogen density field. 
The visual morphology of different models with the same x,n 
values (e.g., Figure|8]) did appear to be noticeably different, 
rather more than one would expect give the tight locus of 
ro — Xm curves. Given the visual impression of the HI slices, 
however, it is understandable that the model with perhaps 
the greatest difference to the others ( the AI > lO^^h^^M© 
model) is also discrepant in terms of ro. For any value of 
Xm, this model has a larger ro than the others. 

The ro differences between models are small, but mea- 
surable, as are the purely gravitational instability based dif- 
ferences (note that the models reach different values of Xm 
at different redshifts, so that the underlying curves will 
be different.) 

We do find a wider variation in the values of the slope 
of ^Hi, for fixed values of Xm- This can be seen from the 
middle panel of Figure 1161 where for example 7 varies from 
1.6 to 0.6 when Xm ~ 0.5. The fiducial model has 7 = 0.72 
when Xm = 0.5. If we look horizontally to find out at what 
Xm value the other models have the same slope, we find 
a range from 1 — Xm ~ 0.34 (for the L/4 model) to 1 — 
Xm = 0.60 (for the M > lO^^h'^M© model). AU curves 
(except for the L/8 model which does not get more than half 
ionized) do follow the same pattern, with the slope getting 
asymptotically fiatter as reionization proceeds, all ending 
up with 7 ~ 0.5. The rapidity with which this asymptotic 
value is reached does vary with the different models, with 
the M > lO^''h~'^M0 model doing this most quickly, as the 
large-scale modulation effects of the bright sources tilt the 
correlation function upwards on large scales. 

The slope of ^hi stays approximately constant as ro 
increases rapidly during the end stages of reionization, as 
can be seen clearly from the bottom panel of Figure 1161 
where ro — 7 curves are plotted. In the region to the left 
of the plot, clustering of HI in all models is dominated by 
the clustering in the underlying density field. The RIS then 
takes over, and again all models follow a rather similar locus 
of £,Hi parameters. We have seen above that the variation 
between the ro — 7 lines is mostly due to the variation of 7 
with Xm in the different models. 

Overall, the change in slope for the different models can 
be qualitively explained by the different length scales of the 
RIS features that occur in each, even at the same stages 
of reionization {xm values). There is a competition between 
small and large scale features that sets the slope of the cor- 
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Figure 17. The growth factor of perturbations in the fiducial 
model. We show results for the gas density field as solid lines 
and the HI density field as dashed lines, with a different line for 
3 different comoving scales. In each case we normalize by the 
amplitude of fluctuations at redshift z = 10.5. The dotted line is 
a curve with g oc e"'*' . 

relation function. The behaviour of ro is more puzzling, and 
we shall return to these questions in Section 6.2 

5.5 The growth factor of perturbations 

We have seen in previous sections that once reionization is 
visually progressing, the amplitude of fluctuations in the HI 
field grows extremely rapidly. It is interesting to compare 
this quantitatively to the growth expected of perturbations 
under gravitational instability. In the latter case, the ampli- 
tude of linear density perturbations will grow at the same 
rate independent of scale, so that ^p(r) will retain the same 
shape, but increase in amplitude by an overall factor [(/(z)]^ 
which can be computed from first order perturbation the- 
ory (e.g., Peebles 1980). At the high redshifts relevant here 
{z ^ 7), the LCDM universe behaves in a manner very close 
to an Einstein de Sitter model so that the linear growth fac- 
tor g{z) (X From Figure [T^ we have seen that ro for the 
matter correlation function is between 0.1 and 0.5 h~^Mpc 
for the models over the range of redshifts when reionization 
occurs, so that linear theory should be accurate over most 
of the range of scales indicated by our resolution and box- 
size tests (Section 5.1). We can also see from Figure [TS] that 
^p(r) does keep the same shape, and increases in amplitude 
slowly, as expected. This can also be inferred from Figure 
1171 where we plot the square root of the ratio of the correla- 
tion function £,p{r) at redshift z to that at redshift z = 10.5, 
as a function of z. This is proportional to the growth factor 
of pertubations between redshifts. There are three lines on 
the plot, corresponding to r = 0.2 h~^Mpc,r = 1 /i~^Mpc 
and r = 7 /i~^Mpc. The lines are close together, indicating 
that the shape of ^p(r) does not change dramatically. 

In Figure [TT] we plot the same quantity for ^///(r), for 
the same r values, for the flducial model of reionization. In 
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Figure 18. The growth factor of perturbations in the HI density 
vs Xm for all 12 models described in §2.3. We show results for 
the 1 /i~^Mpc scale, and have normalized each curve by the 
amplitude of perturbations measured at the last output time to 
be completely neutral. 

this case, we see both that growth of pertubations is much 
more rapid, due to the RIS, than the linear density growth. 
We draw a smooth curve corresponding to g{z) oc e"'*' 
alongside the simulation growth factor for r = 1 /i~^Mpc in 
order to show roughly how extreme the growth in fluctua- 
tions as a function of redshift is during the epoch of reion- 
ization. The growth is faster than exponential over the short 
interval between 2 ~ 10 and z ~ 8, it is roughly an expo- 
nential function of a power law. 

We also see that the different scales exhibit different 
growth rates, with the large-scale fluctuations changing most 
rapidly. On 7 h~^Mpc scales, the fluctuations are first sup- 
pressed, with the square root of ^Hi{r) decreasing by a factor 
~ 5, before rapidly increasing after z = 9.5. This stronger 
behaviour relative to the smaller scales results in the flatten- 
ing of the correlation function. The roughly parallel nature 
of the curves after z ^ 8.5 indicates that at the late stages 
of reionization the flatter power law form of (,Hi{r) has been 
reached, and the amplitude grows similarily on all scales. 

If we now consider the growth of perturbations in the 
different models, a better way to compare them is to look 
at the growth as a function of x,n. In Figure [TSl we show 
this for the r = 1 h~^Mpc scale. As the models become 
more ionized, from left to right, the amplitude of fluctua- 
tions first dips and then rises steeply. The model which dips 
the least on these scales is the M > 10^''h~^MQ model, 
as we expect from looking at Figure 1151 The main growth 
phase of fluctuations starts between x,n = 0.1 Xm = 0.4, 
depending on the model, with an exponential relationship 
between g and Xm- The curves appear to converge towards 
the end of reionization, so that all models exhibit roughly 
the same amount of growth, within ~ 50% from the start of 
reionization until Xm = 0.97. The models which started off 
with less growth in ^hi therefore have steeper dependence 
on Xm- We find that g oc e^'*^™ approximately holds over 
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the range Xm ~ 0.2 — 0.8 for the fiducial model. The steepest 
curve has g cx e"^'^^"" (the L/4 model) and the shallowest has 
g oc e^-^"^'" (for the L > 10^°h-^Mo model). 

If we look at the growth factor from the point of view 
of the parameters varied in each model, we can see that 
decreasing the luminosity (from the fiducial model through 
L/2, L/4 and L/8) monotonically changes how abrupt the 
growth of HI fiuctuations is during the bulk of the reioniza- 
tion process. The model with least luminous sources {L/8) 
exhibits fastest growth with respect to Xm, although of 
course with respect to z, this is not necessarily the case. 
We have examined the growth factor vs z for the difi'erent 
models (not plotted) and find that the growth of fiuctua- 
tions for all models has (at least for the 1 h~^Mpc scale) a 
form roughly consistent wit the g{z) oc e"'*' curve drawn 
on Figure [TT] To zeroth order, the growth of HI fiuctuations 
as a function of z does not seem to be dependent on the 
source physics. 
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5.6 Babul and White model fit 

We have seen that a simple power law fit to the HI correla- 
tion function works well on scales for which the simulation 
has sufficient boxsize and spatial resolution. As reionization 
proceeds, the power law becomes shallower, moving away 
from the slope of the underlying matter correlation func- 
tion. While this paper's focus is on numerical modelling and 
phenomenology of HI clustering rather than analytic mod- 
elling, may neverthless be instructive to consider a different 
fitting function, for ^///(r), derived from very simple model. 
Many analytic models of bubble growth during reionization, 
of varying complexity (e.g., using the excursion set formal- 
ism, Furlanetto et al. 2004, and perturbation theory, Zhang 
et al. 2007.) have been proposed, and most recently, tests 
of the model of Furlanetto, Zhan, and Hernquist have been 
carried out by Zahn et al. (2007) using RT simulations. Here 
we consider the model of Babul & White (1991) originally 
proposed by those authors as a simple description of the 
way galaxy clustering could be modulated by "spheres of 
avoidance" around quasars. 

The model makes the simplifying assumption that the 
sources are Poisson distributed and that they give rise to 
spheres of avoidance of a fixed physical size. The two point 
correlation function of material distribution uniformly in the 
interbubble regions is then (Babul & White 1991): 



1 + Cintcrbubblc(r) 
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(5) 



Here — 4nnbRi/3 is the nominal filling factor of 
spherical bubbles of radius and mean number density nb- 
If we furthermore assume that the distribution of sources is 
independent of the density distribution, modulation of the 
HI by bubbles leads to an HI correlation function given by: 
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(6) 



Of course, the distribution of sources is closely related 
to the density distribution, something that leads to the an- 
tibiasing seen in the early stages of reionization (see Figure 
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Figure 19. A fit of the functional form taken from Babul & 
White (1991) to the autocorrelation function of interbubble ma- 
terial (see Equation |6ll . We show results for the fiducial simula- 
tion, for 8 output times with different values of the mean mass 
weighted ionized H density, Xm ■ The fits were carried out to points 
with r > 0.3 /i-^Mpc. 



I13p . At later times, however, this correlation may become 
less important in governing ^hi, and the simple interbubble 
model may be useful. 

Rather than computing values of the ft and Rt param- 
eters from theory, we use equation [5] to fit ^interbubble mea- 
sured from the simulations, allowing the values of these pa- 
rameters to vary. Before fitting, we first compute ^interbubble 
using the measured value of ^p{r) in the simulation and 
Equation |6l ^interbubble for the fiducial model for various val- 
ues of Xm is shown in Figure [TS] For Xm < 0.4, ^interbubble 
is negative, an unphysical result which occurs because of 
the correlation between sources and p. As reionization pro- 
ceeds, ^interbubble bccomcs dramatically larger in amplitude, 
and the right hand cut off moves to larger scales. The cor- 
relation function is close to flat on smaller scales than the 
break. 

We can interpret this behaviour approximately in the 
context of the Babul and White model fit, which is also 
shown in Figure 1191 We have fit the curves in the same 
manner as was carried out for the power law fits in Sec- 
tion 5.4 (assuming Poisson errors). It is evident that the 
simple model can reproduce the rough shape of ^interbubble, 
with agreement becoming best towards the end stages of 
reionization. It is the modulation of the density correlation 
function by this form for the RIS ^ results in the fiattening 
of ^Hi{r). 

It would be interesting if the best fit value of Rb for a 
given ^interbubble could be used as a measure of the bubble 
size. There are two problems with this, however. The first 
is that in the early stages of reionization when bubbles are 
better defined, they are strongly correlated with sources, so 
that the fit doesn't work. The second problem is that at 
late times, the simulation box size we have used becomes 
comparable to the bubble size and the results will presum- 
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Figure 20. The autocorrelation function of interbubble material 
(see Equation[6ll. We show results for the 12 simulations described 
in §2.3, for the output times when the mean mass weighted ionized 
H density, Xm=0.b. 

ably not be very reliable. There is a range of reionization 
stages for which this might be useful though. For the fits 
in figure [19] the Rb values did vary monotonically with Xm, 
with Rb = 3.2 h-^Mpc for Xm = 0.4, Rb = 7.4 h'^Mpc 
for Xm = 0.5 and Rb ~ 8.1 /i~^Mpc for Xm = 0.7. Above 
Xm = 0.8, Rb saturated at a value of 9.0 /!,~^Mpc, possibly 
indicating the limitations of finite boxsize. 

The other parameter, the filling factor /(, can 
also be examined, and could also be a useful diag- 
nostic. In this case the value steadily increases, with 
fb = 0.06,0.12,0.25,0.46,0.81,1.4,2.7,3.9 for Xm = 
0.4,0.5,0.6,0.7,0.8,0.9,0.97,0.99. We note from the defini- 
tion of fb above that it is a nominal filling factor, and as a 
result fb will become greated than unity as bubbles overlap 
in the late stages of reionization. 

The difference between ^intcrbubbio for the different mod- 
els (all with Xm = 0.5) is shown in Figure 1201 We can see 
that the peak amplitude of the curves varies by a factor 
of ~ 6 between models, but that they all display the same 
shallow curves behaviour with a cutoff. As we would ex- 
pect, the model with the largest obvious bubbles in Figure 
[S] M > lO^''h~^M0 has the largest amplitude and largest 
cutoff scale. A few of the models (for example L/4 and L/8) 
have a negative ^intcrbubbio and are not plotted. For clar- 
ity, we do not show the model fits, but we note that the 
M > 10^° has the largest value of fb, 0.27, and the L/2 
model has the smallest value, 0.06. These values for the bub- 
ble filling factor appear by eye to be at least indicative of 
the trends seen in Figure [S] 

6 SUMMARY AND DISCUSSION 
6.1 Summary 

In this paper we have explored the radiation-induced struc- 
tures (RIS) produced during the epoch of reionization and 



compared them to the results of gravitational instability. 
Our conclusions can be summarized as follows: 

(1) For most of the models we have tested, we find that the 
mean ionized fraction of hydrogen increases exponentially, 
Xm = e^'^"^'' (where Zi is the redshift of full ionization) as 
reionization proceeds. For models adjusted to have emitted 
the same number of source ionizing photons by 2; = 6, there 
is still quite a wide spread in the redshift of reionization, 
with models reaching Xm ~ 0.5 from z — 9.0 to z = 6.4 (the 
last one to reionize is the low erg model.) 

(2) At a fixed Xm, the morphology of the RIS is most 
strongly affected by the lower cutoff in source luminosity, 
which changes the size of bubbles. The mass to light ratio 
of sources also has a substantial effect, but the recombina- 
tion rate and the amplitude of mass fluctuations, a only 
minimally change the appearance of the HI density fleld. 

(3) The HI correlation function, ^hi exhibits a generic 
behaviour for all models tested. S^hi initially becomes 
linearly antibiased with respect to the matter ^p, as the 
high density HI around sources is ionized. The linear bias 
factor reaches a minimum of b ~ 0.5 when Xm ~ 0.1 — 0.2. 
The amplitude of ^h/ then increases rapidly, and (,hi keeps 
a scale invariant shape, but the power law index flattens to 
an asymptotic value of 7 ~ —0.5. 

(4) We flnd that ro, the correlation length of ^hi has the 
essentially the same functional relationship with Xm in all 
but one of the models we have tested. How the power-law 
index 7 varies with Xm on the other hand depends much 
more widely on the different source and physics prescrip- 
tions adopted. 

(5) The growth factor of HI perturbations is seen to change 
much more rapidly than that of gravitionally evolving 
matter perturbations over a redshift range Az ~ 2 — 3 
during which the bulk of reionization occurs. We find 
pertubations on a scale of 1 h~^Mpc to be evolving oc e"'*' 
compared to oc a{t) for gravitational growth. This is valid 
for all models tested, so that the source physics does not 
appear to affect relation. 

(6) During the late stages of reionization, the shape 
evolution of ^hi can be approximately reproduced by 
a simple model due to Babul & White (1991) in which 
ionizing sources are uncorrelated with the density field and 
produce spherical bubbles. Fitting the parameters of this 
model to £^hi therefore form a method for inferring sim- 
ple morphological characteristics from measurements of ^hi- 



6.2 Discussion 

In this paper, we have largely avoided discussing directly 
observational probes of the RIS and the reionization epoch. 
This said, many of our findings can be related closely to the 
possible results of a survey of 21cm brightness, as a function 
of angular position and wavelength. We have concentrated 
mainly on the correlation function of the HI, and so in the 
most likely scenario in which there have been enough early 
X-ray sources to heat up the primordial gas just before reion- 
ization occurs, its temperature is higher than the CMB and 
we are in the emission regime. In this case, the 21cm bright- 
ness temperature is independent of the spin temperature. 
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The brightness temperature T oc a;Hi(l + S), where 5 is the 
overdensity of hydrogen and xm the neutral fraction (see 
e.g., Madau 1997, Di Matteo et al. 2004, Barkana 2007.) 

Looking at the autocorrelation function of the HI, the 
power law behaviour wc sec is striking. For structure growing 
through gravitational instablity, the initial fluctuations are 
scale-invariant over a large range and they evolve in a scale 
free manner (Peebles 1974, Davis & Peebles 1977.) In the 
case of RIS, the onset of significant clustering is extremely 
rapid, and it might be expected that this would lead to fea- 
tures in the correlation function related to for example the 
Stromgren radius of sources dominant at that time. On the 
other hand, cosmic variance would lead to a substantial scat- 
ter in the scale length of features from place to place. Our 
fiducial simulation volume, at 40 comoving /i~^Mpc side 
length also limits our ability to capture the late stages of 
reionization. Wyithe and Loeb (2004) have predicted a co- 
moving radius of bubbles ~ 40 /i~^Mpc at the end of the 
overlap stage. We have seen in §5.6 that the simple Babul 
and White model fit does capture the shape of ^(r) reason- 
ably well even though there is a distinct bubble scale in the 
model. Scale invariance therefore appears relatively easy to 
achieve and detection of departures from it probably requires 
a wide range of scales to be available. 

Purlanetto ct al. (2006) have used the analytic model of 
Purlanetto et al. (2004) to describe how the HII bubble size 
is related to both the bias of galaxies and the underlying 
matter power spectrum. The success of such analytic mod- 
els in comparisons to simulations (e.g., Zahn et al. 2007) has 
opened up the way for their use in analyzing future 21cm 
observations. Our related work, attempting to directly simu- 
late a relatively wide range of models has found that the au- 
tocorrelation function of 21cm emission can be used to infer 
the broad signatures of RIS compared to gravitational struc- 
ture (e.g., non-monotonic growth, flatter asymptotic slope 
7 ~ —0.5) which will help in analysis of the flrst observa- 
tions. 

Among the other analytic models which have been de- 
veloped, the pertubation theory approach of Zhang et al. 
(2007) is difi^erent from many in that it does not make a 
step function bubble approximation to the HI distribution. 
As our simulation approach includes recombinations and the 
contribution of residual HI in the ionized regions, future 
comparisons will be beneficial. For example, Zhang et al. 
(2007) compute the rapid rise in the bias of HI clustering 
as a function of rodshift, finding some qualitatively similar 
results to ours, although they find a scale-independent bias 
on large scales. 

One aspect which we have not covered is the effect of 
redshift distortions on the RIS. We expect the autocorrela- 
tion function to be affected by the usual line of sight ampli- 
fication (Kaiser 1987) on large scales, and small scale sup- 
pression from the velocity dispersion (e.g. Peebles 1980) on 
small scales. The latter effect in particular is likely to be 
strongly modified by the fact that most of the HI around 
bright sources and hence in dense regions is ionized early 
on. These effects will likely be important in the use of 21cm 
emission maps to carry out tests of cosmic geometry (see 
e.g., Nusscr 2005.) It would be simple enough in future work 
to look at the simulations in redshift space, such as would 
be seen with observational data. 

In carrying out our numerical experiments, we have sim- 



ulated a range of models, which we expect to have many of 
the features likely in most scenarios for the reionization of 
the Universe. It was not possible to be completely general, 
however, and it is certainly possible to imagine other inter- 
esting models and tests of the physics that could be included. 
For example, McQuinn et al. (2007) ran as one of their many 
models one in which all cells in the simulation were set to the 
mean density. This had the effect of changing the amount of 
structure in the ionization fronts, and of course drastically 
reducing the recombination rate. 

With our postprocessing RT carried out on the hydro- 
gen distribution only, by keeping track of the ionization state 
but not predicting the temperature wc necessarily simplified 
much of the physics involved in reionization. As our inten- 
tion was to capture the broad differences between RIS and 
gravitational structure, we do not regard this as important. 
However in future work extending our simulation approach 
so that it is directly relevant to upcoming observational 
data, we plan to capture more detail. For example, the code 
SPHRAY (Altay et al. 20007) which represents an extension 
of the present method includes different ionized species and 
explicit temperature evolution. Other work such as Ciardi 
et al. (2006) and McQuirm et al. (2007) specifically include 
the effect of minihalos, sources below the resolution scale 
of the simulation. The latter group finds that they have a 
strong effect on the growth of large bubbles in the late stages 
of reionization. In the present work, we have seen from our 
resolution studies that the autocorrelation function at least 
is not sensitive to increases in mass resolution and small 
scale structure, at least to the accuracy and range of models 
that we have considered. In principle, the high spatial res- 
olution of our gridless ray tracing approach could allow the 
effects of Lyman-limit systems to be modelled, which bo- 
come dominant in limiting the photon mean free path when 
the universe is mostly ionized (e.g., Miralda-Escude et al. 
2000.) 

Future work on realistic models should also include ra- 
diative cooling in the formation of sources and modelling of 
specific sources and spectra. Quasars and miniquasars have 
been investigated by including the formation and growth of 
black holes together with a model for feedback directly into 
cosmological SPH simulations (Di Matteo et al. 2007, Si- 
jacki et al. 2007, Pelupessy et al. 2007). Using such models 
as sources in RT calculations would enable us to investigate 
how the RIS caused by harder sources is different to softer 
sources, in ways which arc not constrained by the simple 
association of source and halo mass as has been carried out 
here. 
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